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Abstract 

Dynamic  characteristics  and  performance  of  a  PEM  fuel  cell  stack  are  crucial  factors  to  ensure  safe,  effective  and  efficient  operation.  In  particular, 
water  and  heat  at  varying  loads  are  important  factors  that  directly  influence  the  stack  performance  and  reliability.  Herein,  we  present  a  new  dynamic 
model  that  considers  temperature  and  two-phase  effects  and  analyze  these  effects  on  the  characteristics  of  a  stack. 

First,  a  model  for  a  two-cell  stack  was  developed  and  the  simulated  results  were  compared  with  experimental  results.  Next,  a  model  for  a  20-cell 
stack  was  constructed  to  investigate  start-up  and  transient  behavior.  Start-up  behavior  under  different  conditions  where  the  amplitudes  and  slopes  of 
a  load  current,  the  temperature  and  flow  rate  of  the  coolant,  and  extra  heating  of  end  plates  were  varied  were  also  analyzed.  The  transient  analyses 
considered  the  dynamics  of  temperature,  oxygen  and  vapor  concentration  in  the  gas  diffusion  media,  liquid  water  saturation,  and  the  variations  of 
water  content  in  the  membranes  at  a  multi-step  load. 

Comparative  studies  revealed  that  the  two-phase  effect  of  water  predominantly  reduces  oxygen  concentration  in  the  catalysts  and  subsequently 
increases  the  activation  over-potential,  while  temperature  gradients  in  the  cells  directly  affect  the  ohmic  over-potential.  The  results  showed  that  the 
heat-up  time  at  start-up  to  achieve  a  given  reference  working  temperature  was  inversely  proportional  to  the  amplitude  of  the  current  density  applied 
and  the  flow  rate  and  temperature  of  the  coolants.  In  addition,  the  asymmetric  profile  of  the  stack  temperature  in  the  stack  was  balanced  when  the 
temperature  of  the  coolant  supplied  was  reheated  and  elevated.  Analyses  of  transient  behaviors  for  a  20-cell  stack  showed  that  strong  temperature 
gradients  formed  in  the  last  four  end  cells,  while  temperature,  oxygen  concentration,  vapor  concentration,  liquid  water  saturation,  and  membrane 
water  content  in  the  rest  of  the  cells  were  uniform. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Polymer  electrolyte  membrane  fuel  cells  (PEMFC)  are  the 
best  candidates  for  an  efficient  power  source  to  supply  our  power 
needs  in  the  future.  This  is  particularly  true  because  of  their 
low  operating  temperature,  relatively  short  start-up  time,  high 
power  density,  and  efficiency.  In  addition,  the  PEM  fuel  cell  pro¬ 
duces  zero  emissions  because  the  chemical  reaction  of  hydrogen 
and  oxygen  ejects  pure  water.  These  advantages  are  superior  to 
those  of  internal  combustion  engines  and  PEMFC  hold  promise 
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for  a  variety  of  applications  in  vehicles,  utilities,  and  portable 
electronics. 

The  life  span  required  for  the  fuel  cell  stack  is  roughly 
3000-5000  h  for  passenger  cars,  up  to  20,000  h  for  commer¬ 
cial  vehicles,  and  up  to  40,000  h  for  stationary  applications.  The 
durability  required  depends  upon  operating  conditions  at  start¬ 
up,  normal  operation,  and  shutdown.  Particularly,  heat  and  water 
management  determine  the  durability  of  the  catalysts.  For  exam¬ 
ple,  improper  operating  conditions  such  as  low  reactant  flows, 
high  and  low  humidity,  or  high  and  low  temperature  can  accel¬ 
erate  Pt  particle  growth  size  and  dissolve  the  particles  from  the 
carbon  support  [1-4].  As  a  result,  the  electro  catalyst  surface 
area  decreases  and  the  overall  performance  drops  because  of 
the  growth  in  platinum  particle  size.  Not  only  a  low  humidity 
but  also  an  elevated  temperature  contributes  to  particle  growth 
[2,3].  It  should  be  emphasized  that  water  and  heat  in  an  oper- 
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Nomenclature 

a  water  activity 

A  area  (m2) 

b  parametric  coefficient  (V) 

c  concentration  of  species  (mol  cm-3) 

C  mas  s  concentration  (kg  m-  3 ) 

Cp  specific  heat  capacity  (J  kg- 1  K- 1 ) 

Dm  diffusion  coefficient  (m2  s-1) 

(An )  effective  diffusivity  (m2  s  - 1 ) 

f^water  water  diffusion  coefficient 

E  open  circuit  potential  (V) 

i  current  density  (A  m-2) 

I  current  (A) 

m  mass  (kg)  or  O2  of  vapor 

M  molecular  mass  (kg  mol- 1 ) 

n  number  of  electrons  (2  or  4) 

rid  electro-osmotic  drag  coefficient 

p  pressure  (atm) 

Q  heat  transfer  ratio  (W) 

Re ie  electrical  resistance  (£2) 

s  liquid  water  saturation  ratio 

AS  change  in  entropy  (J  mol- 1  K- 1 ) 

t  thickness  (cm  or  m),  time  (s) 

T  temperature  (K) 

V  volume  (m3)  or  voltage  (V) 

W  mass  flow  (kg  s-1) 

Greek  letters 
s  porosity 

§  constant  parametric  coefficient 

r)  over-potential  (V) 

X  water  content 

p  mass  density  (kg  m-3) 

(Tmem  membrane  conductivity  ( Q ~ 1  cm- 1 ) 

Subscripts 
an  anode 

ca  cathode 

i  cell  number 

in  inlet 

j  layer 

k  reactant  or  product  gases 

1  anode  or  cathode 

mem  membrane 

sat  saturation 

sou  source 

v  vapor 


ating  stack  are  the  key  factors  that  determine  reliability  and 
durability. 

To  date,  most  of  the  internal  variables  in  an  operating  stack  are 
inaccessible  and  hard  to  measure  because  the  cells  are  composed 
of  thin  layers  and  fully  sealed.  Moreover,  the  interrelated  vari¬ 
ables  in  an  operating  cell  as  well  as  a  stack  make  it  difficult  for 


the  design  and  system  engineers  to  understand  the  mechanism 
and  to  optimally  design  the  stack.  Therefore,  simulations  have 
become  more  accepted  methods  for  studying  the  stack  behavior. 
This  requires  computer  models  that  describe  the  real  behaviors 
of  the  fuel  cell  stack. 

Currently,  models  for  the  PEM  fuel  cell  [5-10]  can  be  catego¬ 
rized  into  multi-dimensional  and  one-dimensional  models.  The 
multi-dimensional  model  solves  governing  equations  by  using 
Computational  Fluid  Dynamics  and  provides  a  high  resolution 
of  flow  characteristics  of  reactants,  fuel,  byproducts,  and  trans¬ 
port  of  charges.  This  is  particularly  advantageous  and  useful  for 
analysis  of  a  single  cell,  but  is  limited  for  a  stack  because  of 
the  enormous  computational  time  and  the  associated  parameters 
and  variables  required  for  the  models.  Moreover,  it  is  impossi¬ 
ble  to  integrate  system  components  into  a  stack  model.  Thus, 
these  models  are  usually  applied  to  investigate  parts  of  complex 
domains  such  as  flow  fields  of  a  single  cell  and  to  predict  the 
performance  of  a  single  cell. 

By  contrast,  the  one-dimensional  model  or  quasi-one- 
dimensional  model  is  a  macroscopic  model  that  describes  layers 
in  a  cell.  These  simplified  models  do  not  provide  detailed  ana¬ 
lytical  mechanisms  of  a  cell,  but  allow  representations  of  the 
dynamic  behavior  of  a  multi-cell  stack  and  analysis  of  the  inter¬ 
actions  between  the  stack  and  system  components.  Therefore, 
these  models  are  preferred  for  investigating  the  dynamics  of  the 
system. 

Many  authors  have  proposed  such  models.  Springer  et  al. 
[10]  proposed  a  model  capable  of  analyzing  water  transport  in 
the  membrane  and  the  performance  of  a  cell  even  though  it  was 
isothermal  and  one-dimensional.  Nguyen  and  White  [11]  devel¬ 
oped  a  model  used  for  analysis  of  the  effects  of  humidification 
and  temperature  on  cell  performance.  Amphlett  et  al.  [12]  pro¬ 
posed  empirical  equations  representing  I-V  characteristics  of 
a  cell  at  different  working  temperatures.  However,  none  of  the 
models  were  able  to  represent  the  dynamics  of  the  physical  prop¬ 
erties  inside  the  cell  at  a  varying  load.  Pukrushpan  et  al.  [13]  pro¬ 
posed  a  model  that  considered  the  dependence  of  the  proton  con¬ 
ductivity  on  the  water  concentration  and  temperature.  However, 
the  water  concentration  of  the  membrane  obtained  from  the  rela¬ 
tive  humidity  (RH)  was  calculated  from  an  average  of  the  anode 
and  cathode  RH.  In  fact,  the  amount  of  water  in  the  membrane  is 
larger  than  that  residing  in  the  anode  and  cathode  side.  Thus,  the 
RH  in  the  anode  and  cathode  varied  rapidly,  while  the  RH  in  the 
membrane  varied  slowly  [14].  In  addition,  the  oxygen  concentra¬ 
tion  in  the  gas  diffusion  layer  (GDL)  on  the  cathode  side  was  con¬ 
tinuously  changing  in  operating  environments  and  significantly 
affected  the  performance  of  the  cell.  Pathapati  et  al.  [15]  added  a 
double  layer  capacitance  that  was  thought  to  represent  the  effects 
of  charges.  However,  both  models  assumed  that  the  behavior  of 
a  stack  was  the  cell  number  times  the  single  cell  performance 
and  that  the  working  temperature  in  a  stack  was  constant. 

A  one-dimensional  thermal  model  developed  by  Khandelwal 
et  al.  [16]  was  used  to  predict  the  temperature  distribution  at  a 
cold  start.  However,  the  model  did  not  provide  transient  analysis 
of  temperature  at  a  varying  load.  Explicit  thermal  dynamic  mod¬ 
els  with  the  coolant  channel  coupled  with  a  ID  cell  model  were 
proposed  by  Amphlett  et  al.  [17],  Wohr  et  al.  [18],  and  Wetton 
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Fig.  1.  Experimental  apparatus. 


et  al.  [19].  In  addition,  Sundaresan  and  Moore  [20]  significantly 
complemented  dynamics  using  simplified  thermodynamic  mod¬ 
els  to  analyze  the  performance  of  the  stack.  The  model  was  based 
on  layers  and  facilitated  analysis  of  start-up  behavior  from  a  sub¬ 
freezing  temperature.  Effects  of  a  varying  load  on  dynamics  of 
temperature  were  also  investigated  by  Shan  and  Choe  [21,22]. 
However,  all  of  the  above  proposed  models  assumed  that  water 
in  the  stack  was  vapor  or  ignored  transient  behavior.  McKay  et 
al.  [23]  developed  a  simple  two-phase  model  for  a  gas  diffu¬ 
sion  layer  and  used  it  to  represent  the  effects  of  liquid  water  and 
vapor  on  the  performance  of  the  cells,  del  Real  et  al.  [24]  pre¬ 
sented  a  two-phase  dynamic  model.  However,  neither  of  these 
models  considered  the  dynamic  variation  of  temperature  and  its 
effects  on  the  performance  of  a  stack  and  water  balance  in  the 
membrane. 

The  model  we  developed  is  effective  and  useful  for  describ¬ 
ing  dynamics  at  varying  loads  and  includes  the  following  major 
improvements  over  other  models:  (1)  temperature  distribution 
in  the  through-plane  of  a  cell  and  its  effects  on  the  character¬ 
istics  of  the  membrane,  catalysts,  gas  diffusion  layers,  and  gas 
flow  channels;  (2)  the  dynamic  water  balance  in  the  membranes; 
(3)  the  two-phase  effects  in  gas  diffusion  layers  and  gas  chan¬ 
nels;  and  (4)  a  coupled  20-cell  stack.  This  analysis  included 
comparisons  between  a  single-phase  and  a  two-phase  model, 
start-up  behavior  of  a  20-cell  stack,  and  the  transient  behavior 
at  a  step  load.  Temperature  distribution  and  the  change  of  over¬ 
potentials  along  with  the  associated  cell  voltages  were  the  major 
parameters  investigated. 

2.  Test  equipment 

The  experimental  apparatus  is  shown  in  Fig.  1 .  The  test  equip¬ 
ment  was  devised  to  test  a  PEM  fuel  cell  system  that  consists  of 
a  two-cell  stack  and  the  balance-of-the-plant  that  includes  an  air 


supply,  a  humidifier,  a  hydrogen  delivery  system,  a  coolant  cir¬ 
cuit,  controls,  and  an  E-load.  A  commercial  software  package, 
Lab  VIEW,  is  used  to  correct  and  process  the  data.  The  stack  was 
constructed  with  two  cells  separated  by  a  thermally  conductive 
plate  in  order  to  minimize  the  potential  influence  of  the  coolants 
on  the  working  temperature.  Other  components  used  in  the  cells 
are  the  same  as  those  in  a  typical  assembly.  The  cell  fabricated 
had  an  active  area  of  140  cm2.  The  maximum  electric  power  of 
the  stack  was  80  W. 

3.  Dynamic  stack  model 

3.1.  Model  description 

The  model  was  developed  on  the  basis  of  layers  in  a  unit 
cell  that  consist  of  a  membrane-electrode-assembly  (MEA),  two 
catalyst  layers,  two  gas  diffusion  layers,  and  two  gas  channels 
sandwiched  between  coolant  channels  as  shown  in  Fig.  2.  The 
input  variables  for  the  model  were  current  load,  mass  flow  rate, 
the  gas  components  fraction,  temperature,  pressure,  and  relative 
humidity  of  reactants  as  well  as  the  temperature  and  flow  rate  of 
coolants  at  the  inlets. 

The  main  assumptions  used  to  develop  the  model  are  as  fol¬ 
lows: 

•  Reactants  are  ideal  gases. 

•  There  is  no  pressure  gradient  between  the  anode  and  cathode 
side  (gas  is  transported  by  diffusion,  not  convection). 

•  There  is  no  gas  pressure  drop  between  the  inlet  and  the  outlet 
of  the  gas  channel. 

•  The  temperature  gradient  is  linear  across  the  layers  in  a  fuel 
cell. 

•  The  thermal  conductivity  for  the  materials  in  a  fuel  cell  is 
constant. 
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Unit  Cell 


(D 

Anode  End  Plate 

© 

Membrane 

© 

Anode  Bus  Plate  (Current  Collector) 

© 

Cathode  Gas  Diffusion  Layer 

Cooling  Channel  Plate 

© 

Cathode  Gas  Channel  Plate 

@ 

Anode  Gas  Channel  Plate 

© 

Cathode  Bus  Plate  (Current  Collector) 

® 

Anode  Gas  Diffusion  Layer 

Cathode  End  Plate 

Fig.  2.  Schematic  simulation  domain  for  a  20-cell  stack. 


•  Anodic  over-potential  is  negligible. 

•  There  is  no  current  density  gradient  across  the  cathode  catalyst 
layer  (the  reaction  is  complete  at  the  cathode  catalyst  layer 
surface). 

•  Latent  heat  during  a  phase  transition  is  not  considered. 

Based  on  these  assumptions,  we  developed  the  models 
described  in  the  following  sections.  Static  behaviors  of  a  cell 
model  are  based  on  voltage  equations  and  the  effects  of  temper¬ 
ature  distribution  in  a  cell,  water  balance  in  the  membrane,  and 
gas  dynamics  in  the  cathode  gas  diffusion  layer  with  a  two-phase 
phenomenon  are  described  below.  A  cell  is  constructed  by  the 
connection  of  individual  models  for  layers  as  shown  in  Fig.  2. 

3.2.  Governing  equations 


Table  1 

Parameters  used  for  models 


Thickness 

(m) 

Thermal 
conductivity 
(W  m-1  K-1) 

Density 

(kgm-3) 

Specific  heat 
(Jkg-'K-1) 

End  plate 

0.028575 

0.228 

1800 

1416 

Coolant  plate 

0.001524 

95 

1780 

935 

Gas  channel  plate 

0.001524 

105 

1820 

935 

Diffusion  media 

0.0004064 

65 

840 

558 

Catalyst  layer 

0.0000035 

0.2 

770 

387 

Membrane 

0.000035 

0.21 

1967 

1100 

Active  area  (cm2) 

140 

The  open  circuit  potential  Ei  is  derived  from  the  energy  bal¬ 
ance  between  chemical  and  electrical  energy  [12]: 


A  model  for  a  single  cell  and  stack  was  developed  and  the  =  1.229  —  0.85  x  10  (7)  —  298.15)  +  4.3085 

composition  of  the  individual  layers  and  the  performance  were 
analyzed  in  detail  [21,22].  Table  1  summarizes  the  parameters 
of  the  stack  components. 


ln((#,H2)  +  -  HPi,o2) 


(2) 


3.2.1.  I-V  characteristic 

The  I-V  characteristic  is  given  by  the  difference  between 
the  open  circuit  voltage  Ei ,  and  the  over-potentials  that  include 
the  activation  over-potential  in  the  catalyst  on  the  cathode  side 
^act,  and  the  ohmic  over-potential  in  the  membrane  r^0hm,  and 
the  concentration  over-potential  ry^on- 

Vi  =  Ei  —  ?7z',act  —  0i,  ohm  —  Oi,  con  (1) 


The  activation  over-potential,  ly^ ct  is  expressed  as  a  func¬ 
tion  of  the  working  temperature  and  the  oxygen  concentration 
[12,28,29]: 

7?z,act  =  §1  +  §2 Ti  +  §3  7)  ln(ci  o2)  +  §47)  ln(7)  (3) 

where  §  represents  constant  parametric  coefficients,  7  is  current 
(A),  and  co2  is  concentration  of  oxygen  in  the  catalyst  interface 
of  the  cathode  (mol  cm-3)  (Table  2). 
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Table  2 

Electrochemical  parameters  [12,28,29] 


Parameters 

Value 

ti 

-0.934 

0.00342 

P 

0.00008 

64 

-0.00019 

b 

0.05 

*max 

1.5 

Hence,  the  internal  energy  source  is  mainly  composed  of  the 
entropy  loss  and  the  chemical  energy  required  for  oxygen  and 
protons  to  overcome  the  barrier  of  the  over-potentials  in  both 
catalyst  layers.  In  addition,  other  heat  sources  are  associated 
with  ohmic  losses  by  transport  of  electrons  and  protons  in  a  cell. 


Gsou  —  l  A 


(  TtAS  .  \ 

( - - ^  Vi  ^  ^ele  ) 


(9) 


The  ohmic  over-potential  r)i, ohm  is  obtained  from  Ohm’s  law 

[10,11]: 


where  the  entropy  change  is  AS  =  0.104  Jmol  1  K  1  for  the 
anode,  and  AS  =  —323.36  J mol-1  K-1  for  the  cathode  [16,27]. 


.  _  1  Wm 

hi, ohm  =  l  ^/,ohm  — 

&i,  mem 


(4) 


where  i  is  the  actual  current  density  (A  cm-2),  and  Rohm ,  Wm 
and  crmem  denote  ohmic  resistance  (£2  cm2),  membrane  thick¬ 
ness  (cm)  and  conductivity  (l£2-1cm-1),  respectively.  The 
membrane  conductivity  is  given  as  a  function  of  the  working 
temperature  [10,11]: 


cr, 


=  cr  exp 


1268 


1  1 
303  _  7 } 
cr  =  0.0051 3 9A^mem  -  0.00326 


(5) 


where  is  the  membrane  water  content. 

The  concentration  over-potentials,  rji, con  are  approximated  by 
the  following  equation  [28,29]: 

W, con  =  -bln  (  1  -  —  j  (6) 

V  /max  J 

where  b  is  a  parametric  coefficient  (V)  that  was  derived  from 
the  single  cell  characteristic  and  assumed  to  be  identical  for  all 
of  the  cells  in  the  stack,  i  and  /max  denote  the  actual  and  the 
maximum  current  density  (A  cm-2). 

Finally,  the  stack  voltage  is  the  sum  of  all  of  the  individual 
cell  voltages. 


Vstack  —  ^  yU 


cell 


(7) 


3.2.3.  Water  balance  in  the  membrane 

Water  content  in  the  membrane  determines  the  proton  con¬ 
ductivity.  The  dynamics  of  the  water  content  is  described  by 
two  effects,  the  electro-osmotic  driving  forces  because  of  the 
electrochemical  potential  difference  between  the  cathode  and 
anode,  and  the  diffusion  caused  by  the  water  concentration  gra¬ 
dient  at  the  two  boundaries.  Considering  the  water  mass  flows 
at  the  boundaries  of  the  membrane  layer,  the  dynamics  of  the 
water  concentration  in  the  membrane  was  described  as  follows 
[10,21]: 


G,H20,mass/47H20 


Pdry,n 


—  0.0126Q5h2c>  ,mass  /47h2o 


m 


z,mem 


3d  mem 
d(C„H2o  ,mass  A  /mem) 
d  t 


(10) 


fk/,ele,mem,an  Ik/,ele,mem,ca  +  W/?diff,mem,an  T  IT/?diff,mem,ca 


where  C is  the  mass  concentration  (kg  m-3),  Mis  the  mole  mass 
(kg  mol-1),  p  is  the  membrane  dry  density  and  A  is  the  fuel  cell 
area  (m2). 

The  electro-osmotic  driving  force  exerted  by  two  different 
electrochemical  potentials  at  the  anode  and  cathode  determines 
the  water  mass  flows  of  W^eie,mem,an  and  WiA ie,mem,ca  at  the 
boundaries  of  the  membrane  layer.  In  addition,  the  diffusion 
caused  by  the  water  concentration  gradient  at  the  two  bound¬ 
aries  makes  up  the  mass  flows  of  IT/,diff,mem,an  and  Wi>diff>mem,Ca- 
Those  relationships  are  described  by  Eqs.  (11)  and  (12)  [9,10]. 


3.2.2.  Energy  balance 

If  a  cell  is  assembled  with  cubical  layers  in  which  thermo¬ 
physical  properties  are  isotropic  and  constant,  the  total  energy 
changes  in  a  controlled  volume  equal  the  sum  of  the  energy 
exchange  at  boundaries  and  internal  energy  resources  accord¬ 
ing  to  the  energy  conservation  equation.  In  fact,  the  energy 
exchanges  at  boundaries  occur  in  three  ways:  (a)  mass  flow  into 
each  volume;  (b)  conduction  across  the  cell;  and  (c)  convection 
occurring  between  bipolar  plates  with  the  coolant,  the  reactants 
and  water.  Thus,  the  thermal-dynamic  behavior  can  be  described 
with  the  following  energy  conservation  equation  [22] : 

E-,Cp,^  =  EwkCPt(r'"  -  tj) 

mass  flow  in 

+  G  convection  +  G conduction  H-  Gsou  (8) 


fk/,ele,mem,l  —  ^watcr  A  n^  jy 

F 


(ID 


lk/,diff,mem,l  —  ^ water  A  7) watery,  1 


(6y,  1  6)', mid) 


(12) 


where  the  electro-osmotic  drag  coefficient  n^fb  the  water  con¬ 
centration  Cy  and  the  water  diffusion  coefficient  D water,/,/  are 
calculated  from  the  empirical  equations  [10]. 


=  0.0029A.?,  +  0.05A./J  -  3.4  x  10"19 


(13) 


^  Pdry,mem . 

Ei,  i  =  — ~ - M,  1, 

IVlme.m 


_  Pdry.mem  .. 

w,mid  —  ..  Ai,mem 

IVlme.m 


T^water ,/,i  —  7)^/jexp  l  2416  l 


303 


1  mem,  i 


(15) 
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Dx,i,  1 


10“6 

1(T6(1  +  2(XiA  -  3)) 
10“6(3  -  1.67(A.;,i  -  3)) 
1.25  x  10“6 


2  ^  7./  ] 

3  >  Ajj  >  2 
4.5  >  >  3 

>  4.5 


(16) 


The  boundary  water  content  Ag  is  a  function  of  water  activity 
a/j,  which  is  calculated  from  the  water  vapor  partial  pressure: 

{  0.043  +  17.81^-j  —  39.85 cYiX  +  3 6a\x  1  >  a\\  >  0 

14  +  1.4(014  —  1)  3>flU>  1 

16.8  >  3 

(17) 


Py,i,  1 

T^sat,/,! 


(18) 


3.2.4.  Gas  diffusion  layer  with  a  two-phase  effect 

The  phase  of  water  in  the  GDL  plays  a  significant  role  in 
transport  of  water  and  reactants.  The  model  for  the  GDL  con¬ 
siders  the  mass  transport  in  a  gas  and  liquid  phase  and  the  phase 
transition  between  liquid  water  and  water  vapor.  Most  dynamic 
models  proposed  for  the  GDL  assumed  that  water  generated  and 
supplied  were  water  vapor,  so  that  the  effects  of  liquid  water  was 
generally  neglected.  However,  when  liquid  water  is  involved, 
diffusion  characteristics  are  different  from  those  in  the  vapor 
state.  The  diffusion  in  the  GDL  can  be  redefined  by  introducing 
an  effective  diffusivity  that  describes  the  diffusion  behavior  of 
vapor  and  liquid  water  in  a  capillary  [23,25] .  Porosity  is  reflected 
in  the  equation,  while  no  tortuosity  was  considered. 


{Dm(k))  —  Dms 


£  —  0.11 


0.785 


£  = 


y  pore 


1-0.11 

^liquid  (k) 


(i  -  s(k)? 


Vgdl’ 
m  =  O2,  vapor 


s{k)  = 


Vr 


pore 


(19) 


where  (Dm)  is  the  effective  diffusivity  (m2  s-1),  Dm  is  the  dif¬ 
fusion  coefficient  (m2  s_1)  at  a  single  phase,  £  is  the  porosity  of 
the  diffusion  layer,  5  is  the  liquid  water  saturation  ratio,  Vp0re 
is  the  pore  volume  of  the  GDL  (m3),  k  is  the  number  of  each 
domain  in  the  GDL,  Vgdl  is  the  total  volume  of  the  GDL  (m3), 
and  Viiquid  is  the  volume  of  the  liquid  water  (m3). 

The  simulation  domain  for  half  a  cell  is  depicted  in  Fig.  3, 
where  the  GDL  was  divided  into  three  sub-domains  in  order 
to  calculate  gradients  of  species.  The  points  of  0  and  4  are  the 
boundary  conditions,  while  1,  2  and  3  are  the  middle  points  of 
the  sub-domains. 

The  water  vapor  produced,  Av,o,  and  the  oxygen  consumed, 
No2,0>  were  expressed  as  a  function  of  a  current.  The  total  molar 
flux  of  the  water  vapor  at  the  boundary  was  given  as  follows: 


Nv’°  =  YF’  n°^°  =  4 J 

Ny, catalyst  =  H”  Ny  ,mem 

where  /Vv;mem  is  defined  in  the  Eq.  (10). 


(20) 


Cathode  GDL 


Fig.  3.  Schematic  simulation  domain  for  the  GDL  with  a  two-phase  phe¬ 
nomenon. 


The  molar  flux  is  expressed  as  a  function  of  the  effective  diffu¬ 
sivity  and  the  gas  concentration  gradient,  where  zca(&)  represents 
the  ratio  of  the  species  molar  flux: 


Zca(*)  = 


Noffk)  = 


Ny(k)  = 


)  Nv,o/No2,o  for  *  =  1 
[  Nffk  -  l)/No2(k  -  1)  for  k  =  2,  3 

-(Po2(k))  cp2(k  +  1)  -  cp2(k ) 

1  -  X02(k)(l  +  ZcaW)  ytG DL 

—  (Dy(k))  Cy(k  +  1)  -  Cy(k) 

1  -  *v(fc)(l  +  1/ZcaW)  ytG  DL 


(21) 


N  is  the  molar  flux  (mols-1  m-2),  v  is  the  molar  ratio  of  the 
species,  and  y  is  the  coefficient  dependent  upon  the  number  of 
the  sub-domain. 

The  time  derivatives  of  the  gas  concentrations  are  a  function 
of  the  ratio  of  the  species  molar  flux  [23]: 


dcp2  _  Np2(k)  -  Nq2 (k  -  1) 
d  t  ytGDL 


dc,  Ny(k)  —  Ny(k  —  1) 

—  (k)  = - 

d  t  ytc  dl 


+  ^evap(^) 


(22) 


Hence,  the  evaporation  rate  is  determined  by  the  difference 
between  the  saturation  and  vapor  pressure  [23,25]: 


^evap(^)  —  Y 


Pv,  sat(T^)  Pv(k) 

RT 


(23) 


where  y  is  the  volumetric  condensation  coefficient  (s-1),  p  is 
pressure  (Pa)  and  R  is  the  ideal  gas  constant  (J  mol-1  K-1). 

According  to  the  mass  balance,  the  rate  of  the  liquid  water 
volume  at  each  of  the  points  was  obtained  as  follows: 


d  Viiquid 

d  t 


MwsVGDLRowap(k)  WiiquidC^)  ^  j 

Pliquid 

(k)  =  ^  —  MyS  VGDL^evap(^)  +  WiiquidfA:  —  1)  —  VVliquidf/c) 


(24) 


Pi  i  quid 


k  =  2,3 


The  mass  flow  of  the  liquid  water  in  the  Eq.  (24),  Illiquid, 
was  derived  from  the  pressure  gradient  driven  flow  in  a 
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Fig.  4.  Comparison  of  the  simulated  and  experimental  stack  and  cell  voltages  at  a  two-cell  stack. 


capillary  [23,25]: 

^Pliquid^^rl 


^liquid  (JO  — 


S(k)  = 


d  pc 


dS 


/4  liquid 
S(k)  -  ^jm 
1  —  sim 

0,  0  <  s(k)  <  Sh 


S(k  +  1)  -  S(k) 


ytG  DL 
•Sim  <  s(k)  5  1 


(25) 


where  p  is  the  liquid  water  density  (kgm-3),  K  is  the  absolute 
permeability  (m2),  Kv\  is  the  relative  permeability  of  liquid  water 
that  was  assumed  to  be  equal  to  S'3,  /x  is  the  viscosity  of  water 
(kgm-1  s-1),  \dpc/dS\  is  the  slope  of  capillary  pressure  (Pa),  S 
is  the  reduced  liquid  water  saturation  and  S{m  is  the  immobile 
saturation. 

Table  3  summarizes  the  parameters  of  the  GDL  used  for  the 
simulation.  The  range  of  the  simulated  liquid  water  saturation 
shown  in  Fig.  10  was  0. 12-0. 18,  which  is  comparable  with  those 
in  Ref.  [25]. 


4.  Results  and  discussion 

Simulations  to  analyze  dynamic  performance  of  a  stack  were 
performed.  First,  I-V  characteristics  of  the  stack  as  well  as  the 
individual  cell  for  a  two-cell  stack  were  compared  with  the 
experimental  data.  Second,  the  effects  of  a  single-phase  and  two- 
phase  of  water  on  the  performance  were  evaluated  in  a  two-cell 
stack.  Finally,  a  model  for  a  20-cell  stack  was  constructed  and 
used  to  analyze  dynamic  characteristics  at  a  start-up  and  at  a 
varying  load. 


Table  3 

Parameters  for  the  GDL  [25] 


Parameters 

Value 

£ 

0.5 

Do2 

3.03  x  10“5  (m2  s”1) 

/^vapour 

3.45  x  10“5  (m2  s”1) 

K 

2.55  x  10“ 13  (m2) 

\dpc/dS\ 

30.321  (Pa) 

sim 

0.1 

4.1.  Comparison  between  simulations  and  experiments  for 
a  two-cell  stack 

Using  our  model,  we  simulated  a  two-cell  stack.  The  I-V 
characteristics  calculated  for  the  stack  as  well  as  the  individual 
cells  were  compared  with  the  experimental  data  collected  at  a 
test  station  installed  at  Auburn  University.  Operating  conditions 
were  as  follows  where  the  stoichiometric  ratio  for  the  anode 
was  1.2  and  for  the  cathode  was  3.0.  The  gas  pressure  and  sup¬ 
ply  gas  temperature  on  both  sides  were  1.0 bar  and  333.15  K, 
respectively.  The  relative  humidity  for  the  anode  gas  was  0% 
and  for  the  cathode  gas  was  100%.  The  average  temperature 
of  the  coolants  at  the  exits  of  cells  1  and  2  were  maintained  at 
333.15  K  by  controlling  the  coolant  temperature. 

Fig.  4  shows  simulated  and  experimental  results  of  the  voltage 
across  the  stack  (a)  and  the  cells  (b).  The  simulated  and  exper¬ 
imental  stack  voltages  matched  well,  but  the  deviation  of  the 
simulated  two  cell  voltages  was  not  as  large  as  that  measured 
experimentally.  Differences  in  the  two  cells  could  be  caused 
by  non-uniform  characteristics  of  the  individual  cells  such  as 
reactant  distributions  and  the  properties  of  the  GDL  and  the 
membrane.  Details  describing  the  effects  of  the  various  factors 
on  the  cell  voltages  are  described  in  Park  and  Choe  [26] . 

4.2.  Comparison  of  a  single-phase  and  a  two-phase  model 
for  a  two-cell  stack 

The  structure  of  the  two-cell  stack  was  the  same  as  that 
depicted  in  Fig.  2  except  for  the  addition  of  a  separator  between 
the  two  coolant  channels  that  was  used  to  minimize  the  coupled 
cooling  effects  of  adjacent  cells  on  performance.  Operating  con¬ 
ditions  were  as  follows:  the  stoichiometric  ratio  for  the  anode 
was  1.2  and  that  for  the  cathode  was  2.0.  The  gas  pressure,  sup¬ 
ply  gas  temperature,  and  relative  humidity  on  both  sides  were 
1.0  bar,  333.15  K,  and  100%,  respectively.  The  average  temper¬ 
ature  of  coolants  at  the  exit  of  cells  1  and  2  for  the  single-phase 
model  was  maintained  at  333.15  K  by  controlling  the  coolant 
inlet  temperature,  and  the  coolant  flow  rate  was  constant  at 
0.003 kgs-1.  The  two-phase  model  used  the  same  coolant  tem¬ 
perature  and  flow  rate. 

Static  and  dynamic  characteristics  of  a  two-cell  stack  are 
illustrated  in  Fig.  5.  Fig.  5(a)  shows  two  major  effects  in  that 
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Fig.  5.  I-V  curve  and  dynamic  response  of  fuel  cells  during  step  changes  in  current  density,  (a)  Calculated  polarization  curves  of  single-phase  and  two-phase 
models  for  a  two-cell  stack,  (b)  Variation  of  cell  voltage  during  step  changes  in  current  density,  (c)  variation  of  activation  and  ohmic  over-potentials,  (d)  temperature 
distribution  in  a  two-cell  stack. 
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Fig.  6.  Temperature  distribution  in  the  layers  of  a  20-cell  stack  for  start-up  in  different  cases. 
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Fig.  7.  The  dynamic  response  of  the  cell  voltage  during  step  change  in  current  density  (a),  (b)  coolant  flow  rate  and  coolant  temperature,  (c  and  d)  temperature 
distribution  through  the  layer  during  step  changes  in  the  current  density  (0.4,  0.8  and  1.2  A  cm-2). 


the  voltage  of  the  two  models  were  different  from  each  other 
and  that  the  deviation  of  the  voltages  between  cells  1  and  2  was 
larger  as  the  load  current  increased.  The  difference  in  the  cell 
voltages  between  a  single-and  two-phase  models  is  affected  by 
liquid  water  that  is  otherwise  neglected.  As  the  density  of  the 
load  current  is  higher  than  0.4  A  cm-2,  the  influence  of  liquid 
water  on  the  performance  of  the  cells  increased  and  the  differ¬ 
ence  between  the  cells  was  larger.  As  a  result,  the  voltage  drop 
at  the  two-phase  model  is  larger  than  that  at  the  single-phase 
model.  In  addition,  the  high  voltage  in  cell  2  is  caused  by  differ¬ 
ent  working  temperatures  that  depend  upon  the  location  of  the 


cells  and  the  coolant  channel.  However,  it  should  be  noted  that 
the  liquid  water  did  not  significantly  affect  the  limiting  current 
density  in  comparison  to  results  of  other  analyses  because  the 
gradient  of  water  concentration,  the  pressure  and  temperature 
drop  along  the  flow  channel  have  not  been  fully  considered  in 
the  modeling  [6,30,31]. 

Fig.  5(b)  shows  the  dynamic  response  of  a  current  load  pro¬ 
file  on  the  cell  voltages.  Three  amplitudes  of  the  current  density 
(0.4, 0.8  and  1.2  A  cm-2)  were  applied  to  observe  the  responses. 
Two  major  phenomena  were  observed.  First,  the  voltages  of 
the  cells  with  a  two-phase  model  were  lower  than  those  with 


Cell  number 


Fig.  8.  The  dynamic  response  of  the  oxygen  concentration  during  step  changes  in  current  density  (a),  (b)  oxygen  concentration  throughout  the  stack. 
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Fig.  9.  The  dynamic  response  of  the  vapor  concentration  during  step  changes  in  current  density  (a),  (b)  vapor  concentration  throughout  the  stack. 
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Fig.  10.  The  dynamic  step  response  of  liquid  water  saturation  at  a  nearby  cathode  catalyst  layer  during  step  changes  of  current  density  (a),  (b)  liquid  water  saturation 
distribution  throughout  the  stack. 


a  single-phase  model.  Second,  the  higher  the  current  density 
the  larger  the  voltage  drop  caused  by  the  two-phase  effect. 
Consideration  of  the  presence  of  liquid  water  in  the  gas  dif¬ 
fusion  layer  changes  the  activation  and  ohmic  over-potential 
shown  in  Fig.  5(c).  The  activation  over-potential  is  influenced 
by  the  liquid  water  and  independent  of  the  location  of  the 
cells.  As  the  liquid  water  in  the  gas  diffusion  layer  increased, 
the  diffusivity  of  oxygen  was  reduced  and  its  concentration 


was  reduced.  Consequently,  the  activation  over-potentials  of 
the  cells  became  higher  and  the  voltages  of  the  cells  became 
lower. 

On  the  other  hand,  ohmic  over-potentials  are  dependent  on 
the  location  of  the  cells  and  are  asymmetrical,  and  also  affected 
by  temperature  in  the  cell.  Fig.  5(d)  shows  a  temperature  distri¬ 
bution  in  the  stack  at  the  given  profile  of  current  density.  When 
liquid  water  was  considered,  the  over-potentials  in  the  cells  were 
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Fig.  11.  The  dynamic  response  of  the  membrane  water  content  during  step  changes  in  current  density  (a),  (b)  membrane  water  content  throughout  the  stack. 
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higher  than  those  in  a  single-phase  model  and  subsequently 
the  associated  heat  produced  was  more  than  that  observed  in 
a  single-phase  model.  As  a  result,  the  working  temperature  in 
the  cells  became  higher.  In  addition,  the  geometrical  distribution 
of  the  heat  sources  in  the  stack  was  asymmetrical.  Accordingly, 
heat  conductivity  from  the  heat  sources  to  the  coolants  was  not 
identical,  which  caused  an  asymmetrical  temperature  profile  in 
that  the  temperature  in  cell  1  was  higher  than  that  in  cell  2.  The 
ohmic  over-potential  was  affected  by  the  temperature  and  water 
content  in  the  membrane,  and  the  temperature  in  the  membrane 
of  cell  2  was  higher  than  that  of  cell  1 .  In  addition,  the  cathode  gas 
channel  temperature  of  cell  2  was  lower  than  that  of  cell  1,  and 
the  saturation  vapor  pressure  was  also  lower,  which  augments  the 
RH  in  cell  2  and  water  content.  Consequently,  the  ohmic  over¬ 
potential  in  cell  2  became  lower  than  that  in  cell  1  and  determined 
the  deviation  of  the  amplitude  of  the  two  cell  voltages. 

4.3.  Analysis  of  a  20-cell  stack 

Start-up,  normal  operation,  and  shutdown  are  processes 
needed  for  operation  of  a  stack.  We  investigated  the  behavior  of 
a  20-cell  stack  at  a  start-up  and  under  normal  operations  using 
the  proposed  model. 

4.3.1.  Start-up 

Dynamic  behavior  of  a  start-up  is  particularly  important  to 
ensure  a  short  start-up  time  along  with  a  high  efficiency.  In  gen¬ 
eral,  the  optimal  operating  temperature  of  a  PEM  fuel  cell  is 
from  333.15  to  353.15  K.  This  facilitates  relatively  high  chemi¬ 
cal  reactions  and  allows  the  water  produced  to  be  easily  removed, 
while  maintaining  the  membrane  at  a  proper  operating  temper¬ 
ature. 

A  cold  stack  after  a  shutdown  can  be  deployed  in  different 
environments,  as  example,  under  a  freezing  temperature  where 
moisture  residing  inside  the  stack  should  have  formed  ices  or, 
alternatively,  above  the  freezing  temperature.  Start-up  behavior 
of  a  frozen  stack  used  for  mobile  applications  has  been  theoreti¬ 
cally  and  experimentally  investigated  [16,20,32-34].  However, 
the  start-up  behavior  of  a  stack  above  the  freezing  temperature 
has  not  been  fully  investigated.  In  fact,  a  low  working  temper¬ 
ature  of  a  stack  decreases  its  performance  because  of  increased 
kinetic  losses,  and  ohmic  and  reactant  transport  losses  caused  by 
a  high  rate  of  condensation.  Thus,  it  is  necessary  to  elevate  the 
working  temperature  of  a  stack  as  quickly  as  possible  to  meet 
the  demands  of  the  required  load  power. 

Different  start-up  strategies  were  compared  in  the  studies 
described  below.  The  optimal  working  temperature  of  the  stack 
was  assumed  to  be  333.15  K.  Operating  conditions  were  as  fol¬ 
lows:  the  stoichiometric  ratio  for  the  anode  was  1.2  and  that 
for  the  cathode  was  2.0.  The  pressure,  temperature  and  rela¬ 
tive  humidity  of  hydrogen  on  the  anode  side  were  set  to  1 .2  bar, 

298.15  K  and  0%,  respectively,  while  those  on  the  cathode  side 
were  1.0 bar,  338.15  K  and  100%.  The  unconsumed  hydrogen 
was  recirculated  to  the  stack.  Table  4  summarizes  the  heat-up 
time  of  a  stack  using  different  start-up  strategies.  The  initial 
temperature  of  the  stack  was  assumed  to  be  298.15  K. 


In  cases  A,  B  and  C,  we  assumed  that  current  density  was 
applied  as  a  step.  The  heat-up  time  became  shorter  as  the  ampli¬ 
tude  of  the  current  applied  became  higher.  However,  selection 
of  a  current  density  should  be  carefully  determined  after  consid¬ 
eration  of  the  maximum  allowed  temperature  and  the  associated 
thermal  energy  for  the  membrane  and  the  catalyst.  In  cases 
D,  E  and  F,  the  current  density  was  increased  with  a  different 
slope,  from  0.2  to  1.0  A  cm-2,  where  no  coolant  was  supplied. 
When  the  current  was  increased  within  30  s,  the  stack  had  not 
achieved  the  set  working  temperature.  It  took  a  minimum  of 
40  s  for  the  stack  to  achieve  the  working  temperature.  Case 
G,  assumed  that  the  initial  temperature  of  the  end  plate  was 

333.15  K.  As  compared  with  case  B,  the  heat-up  time  in  case 
G  was  not  reduced  by  heating  the  end  plates.  Cases  H  and 
I  assumed  that  the  initial  temperature  of  the  end  plate  was 

298.15  K  and  that  of  the  coolant  was  298.15  K  with  a  con¬ 
stant  flow  rate  of  0.0005  and  0.003 kgs-1  that  excluded  any 
heat  exchanger.  Compared  with  case  B,  the  heat-up  time  was 
increased  because  of  the  heat  transfer  from  the  stack  to  the 
coolant  when  the  temperature  of  the  coolant  was  the  same  as 
that  of  the  stack  at  the  initial  state.  Cases  J  and  K  assumed 
different  coolant  flow  rates  of  0.0005  and  0.003  kg  s-1,  respec¬ 
tively,  where  the  temperature  of  the  coolants  was  controlled  at 

333.15  K.  Case  J  allowed  for  a  15.7  s  reduction  of  the  start-up 
time  relative  to  case  G.  In  addition,  in  case  K  the  heat-up  time 
was  reduced  by  35.8  s.  Cases  L  and  M  showed  start-up  times 
similar  to  those  of  cases  J  and  K  regardless  of  the  temperature 
of  the  end  plates. 

Temperature  distributions  in  the  20-cell  stack  under  various 
start-up  conditions  are  summarized  in  Table  4  and  shown  in 
Fig.  6(a-d).  The  graphs  for  the  temperature  profiles  were  cap¬ 
tured  after  the  heat-up  time.  Ohmic  losses  by  the  membrane 
resistances  and  the  heat  released  by  the  chemical  reaction  in 
the  catalysts  caused  the  highest  peaks  among  others  with  an 
asymmetrical  temperature  distribution  throughout  the  cell. 

When  the  amplitude  of  applied  current  was  0.2,  0.4, 
0.8  A  cm-2  (cases  A,  B  and  C),  the  temperature  difference  (A T) 
between  the  layer  with  a  maximum  value  (333.15  K)  and  the 
catalyst  temperature  of  cell  1  became  5.26,  7.42,  and  8.76  K, 
respectively.  The  deviation  (AT)  shown  in  Fig.  6(a)  was  larger 
with  increasing  current.  Because  of  the  high  current  density,  a 
large  amount  of  heat  was  produced  in  a  short  time  and  the  tem¬ 
perature  rose  rapidly  through  the  cells,  particularly  in  the  middle 
cell. 

When  the  coolant  inlet  temperature  was  333.15  K  (cases  F 
and  M),  the  temperature  difference  AT  was  reduced  5.25  and 
2.56  K,  respectively.  The  asymmetric  stack  temperature  distri¬ 
bution  was  minimized  when  supplying  coolants  with  elevated 
temperature. 

4.4.  Transient  response 

In  the  following  section,  transient  behaviors  of  a  20-cell 
stack  were  investigated.  Operating  conditions  were  the  same  as 
those  in  Section  4.3.  The  coolant  inlet  temperature  was  kept  at 

333.15  K  as  shown  in  Fig.  7(b).  The  average  temperature  of  the 
coolant  channel  outlet  temperature  was  controlled  at  343.15  K 
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Table  4 

Heat-up  time  at  various  operating  conditions 


Case 

Constant  current 
density  (A  cm-2) 

Time  for  increasing 
current  density 

0.2-1.0  Acm-2  (s) 

End  plate  initial 
temperature  (K) 

Coolant  flow  rate 

for  each  coolant 
channel  (kg  s-1) 

Coolant  inlet 
temperature  (K) 

Temperature 
difference  (AT) 

Heat-up 
time  (s) 

A 

0.2 

- 

298.15 

No-coolant 

- 

5.26 

125.9 

B 

0.4 

- 

298.15 

No-coolant 

- 

7.42 

54.0 

C 

0.8 

- 

298.15 

No-coolant 

- 

8.76 

21.3 

D 

- 

30 

298.15 

No-coolant 

- 

- 

- 

E 

- 

60 

298.15 

No-coolant 

- 

8.41 

42.9 

F 

90 

298.15 

No-coolant 

- 

8.14 

50.6 

G 

0.4 

- 

333.15 

No-coolant 

- 

6.58 

53.9 

H 

0.4 

- 

333.15 

0.0005 

298.15 

(circulation) 

5.44 

55.2 

I 

0.4 

- 

333.15 

0.003 

298.15 

(circulation) 

2.05 

56.7 

J 

0.4 

- 

333.15 

0.0005 

333.15 

4.69 

38.2 

K 

0.4 

- 

333.15 

0.003 

333.15 

2.35 

18.1 

L 

0.4 

- 

298.15 

0.0005 

333.15 

5.25 

38.2 

M 

0.4 

- 

298.15 

0.003 

333.15 

2.56 

18.1 

by  the  coolant  inlet  flow  rate  as  the  current  density  varied.  The 
coolant  flow  rate  refers  to  that  rate  of  coolant  supplied  to  a 
coolant  channel. 

When  a  step  current  was  applied  as  a  load,  the  voltage  of 
the  cells  varied.  As  seen  in  Fig.  7(a),  no  significant  difference 
between  the  individual  cells  was  observed.  In  the  operation 
of  a  real  stack,  there  are  differences  in  the  I-V  characteris¬ 
tic  of  the  individual  cells  influenced  by  inherent  properties  of 
the  components  as  well  as  the  design  parameters  of  the  flow 
patterns  [31]  that  are  not  considered  in  the  proposed  model. 
Fig.  7(c)  shows  the  stack  temperature  distribution  at  0.4  A  cm-2 
(200  s),  0.8  A  cm“2  (235  s),  and  1.2  A  cm-2  (270  s).  We  find 
that  the  higher  the  current  density,  the  higher  the  amplitude  of 
the  peak  temperature.  The  maximum  temperature  was  1.5-2  K 
higher  than  the  average  temperature  of  the  coolant  outlet  at  a 
current  density  of  1.2  A  cm-2.  In  addition,  the  temperature  in 
the  catalyst  of  cell  1  was  lower  than  that  in  cell  20  because 
of  high  heat  conductivity  to  the  coolant  channels.  The  four 
cells  from  the  anode  and  cathode  end  plate  showed  a  temper¬ 
ature  gradient,  while  the  rest  of  the  cells  exhibited  relatively 
uniform  temperature  distributions.  When  the  current  density 
increased,  the  temperature  of  the  anode  coolant  channel  in 
cell  1  was  reversed  because  of  increased  coolant  flow  rates 
and  a  high  gradient.  However,  the  difference  of  1.5-2  K  in  the 
cathode  catalysts  did  not  affect  the  cell  voltages.  Stack  and 
system  design  should  be  optimized  with  caution  to  ensure  a  uni¬ 
form  temperature  distribution  in  the  stack.  Depending  upon  the 
operating  conditions,  the  stack  can  be  locally  clogged  by  con¬ 
densed  liquid  water  when  the  temperature  in  the  region  suddenly 
drops. 

Fig.  7(d)  shows  the  temperature  distribution  of  the  stack  when 
the  coolant  flow  rate  at  the  anode  side  of  cell  1  and  cathode  side  of 
cell  20  was  reduced  by  1/10  of  the  operating  condition  shown  in 
Fig.  7(c).  The  reverse  effect  in  the  anode  coolant  channel  was  not 
observed  and  the  stack  temperature  was  uniformly  distributed 
in  comparison  to  Fig.  7(c).  Different  coolant  flow  rates  at  the 


endplates  can  be  implemented  by  proper  design  of  the  number 
and  dimensions  of  coolant  channels. 

Fig.  8(a)  shows  the  dynamic  response  of  the  oxygen  con¬ 
centration  at  a  nearby  cathode  catalyst  layer  under  the  same 
condition  as  shown  in  Fig.  7(a-c).  In  general,  more  oxygen 
was  consumed  as  the  load  current  increased  and  the  concen¬ 
tration  became  lower.  In  addition,  the  temperature  influenced 
the  concentrations  of  oxygen  in  the  cell.  Fig.  8(b)  shows  the 
oxygen  concentration  at  a  nearby  cathode  catalyst  layer  of 
the  cells.  The  concentration  in  the  end  cells  was  higher  than 
that  in  the  other  cells  because  of  the  temperature  drop  in  the 
cathode  gas  channel  that  decreased  the  saturation  pressure  and 
increased  the  condensation.  Subsequently,  the  partial  pressure 
of  water  vapor  decreased  and  the  partial  pressure  of  oxygen 
increased. 

Fig.  9(a)  shows  the  dynamic  response  of  the  vapor  concentra¬ 
tion  at  a  nearby  cathode  catalyst  layer.  In  general,  more  vapors 
were  generated  as  the  load  current  increased  and  the  oxygen 
concentration  became  lower.  Fig.  9(b)  shows  the  vapor  con¬ 
centration  at  a  nearby  cathode  catalyst  layer  of  the  cells.  The 
concentration  in  the  end  cells  was  lower  than  that  in  the  other 
cells  because  of  the  influence  of  the  temperature  distributions. 

Fig.  10(a)  shows  the  dynamic  response  of  liquid  water  at 
a  nearby  cathode  catalyst  layer  to  step  changes  in  the  current 
density.  The  liquid  water  saturation  on  the  cathode  catalyst 
layer  followed  the  load  profile.  As  the  load  current  increased, 
more  water  was  produced  in  the  catalyst.  Fig.  10(b)  shows 
the  liquid  water  saturation  at  a  nearby  cathode  catalyst  layer 
of  the  cells.  The  liquid  water  saturation  in  the  end  cells  was 
higher  than  that  in  the  other  cells  because  of  the  temperature 
distributions. 

Fig.  1 1(a)  shows  the  dynamic  step  response  of  a  current  den¬ 
sity  on  the  membrane  water  content,  while  Fig.  11(b)  shows 
the  membrane  water  content  of  the  cells.  As  the  load  current 
increased,  the  membrane  water  content  decreased  because  of 
the  high  water  uptakes  from  the  anode  to  the  cathode  side. 


672 


S.-K.  Park,  S.-Y.  Choe  /  Journal  of  Power  Sources  179  (2008)  660-672 


5.  Conclusions 

We  focused  on  the  development  of  a  dynamic  model  for  a 
20-cell  stack  that  considered  temperature  and  two-phase  effects 
in  the  GDL. 

The  major  contributions  of  this  work  are  summarized  as  fol¬ 
lows: 

•  Voltage  differences  between  a  single-phase  model  and  two- 
phase  model  were  predominantly  affected  by  activation 
over-potentials  caused  by  changes  of  oxygen  concentration, 
while  the  voltage  differences  between  cells  1  and  2  were 
caused  by  temperature  gradients  in  each  layer  and  subse¬ 
quent  ohmic  over-potentials.  However,  the  liquid  water  did 
not  significantly  affect  the  limiting  current  density  when  com¬ 
pared  with  the  results  of  other  analyses  because  the  gradient 
of  water  concentration,  the  pressure  and  temperature  drop 
along  the  flow  channel  were  not  fully  considered  in  the  present 
modeling. 

•  Various  strategies  for  a  start-up  of  a  20-cell  fuel  cell  stack  were 
analyzed.  In  general,  the  higher  the  current  density  and  the 
more  the  coolant  flow  rate  and  temperature  were  increased, 
the  shorter  the  heat-up  time.  However,  the  initial  temperature 
of  the  end  plate  did  not  play  a  significant  role  in  reducing  the 
heat-up  time.  The  asymmetric  distribution  of  the  voltage  was 
minimized  by  supplying  the  coolant,  and  maximized  as  the 
temperature  of  the  coolant  became  higher. 

•  Upon  operation  of  a  20-cell  stack,  the  four  cells  from  each 
end  showed  a  temperature  gradient,  while  the  rest  of  the 
cells  maintained  a  relatively  uniform  temperature  distribu¬ 
tion,  oxygen  concentration,  vapor  concentration,  liquid  water 
saturation,  and  membrane  water  content. 

•  Asymmetric  distribution  of  temperature  was  balanced  by 
reduction  of  the  coolant  flow  rate  at  both  end  cells,  which 
provided  a  better  voltage  distribution  at  dynamic  loads. 
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